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urally arise in the context of a wide range of nonlinear models in different numbers of spatial 
dimensions. We present a numerical method, which we call the adiabatic damping method, designed 
to study such configurations in small lattices. Using 3-dimensional oscillons in models as an ex- 
ample, we show that the method accurately (to a part in 10^ or better) reproduces results obtained 
with static or dynamically expanding lattices, dramatically cutting down in integration time. We 
further present new results for 2-dimensional oscillons, whose lifetimes would be prohibitively long 
to study with conventional methods. 
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I. INTRODUCTION 

The study of long-lived, localized coherent configura- 
tions is of great interest to many areas of physics and 
engineering. In general, these structures arise within the 
context of numerical studies of effective nonlinear field 
models, which may either describe behavior already ob- 
served in the laboratory or conjectured to exist in phe- 
nomena yet to be observed. In contrast to the usual 
solitonic behavior, which is marked by localized time- 
independent configurations, these are time-dependent 
configurations, which nevertheless persist for very long 
times. It is reasonable to suppose that the long life- 
times are due to energy exchange promoted by nonlinear 
coupling between different modes which efficiently sup- 
presses the radiation of energy away from the configu- 
ration. One may collectively call these persistent coher- 
ent field configurations, (PCFCs) in order to distinguish 
them from the usual static coherent configurations which 
characterize solitonic behavior. 

Perhaps the most well-known PCFCs are the 1- 
dimensional breathers which appear during low velocity 
kink-antikink interactions in sine-Gordon and 4>'^ mod- 
els [0. These are bound states characterized by nonlin- 
ear oscillations about the energy minimum (or vacuum) 
of the model, typically one of the minima in a degen- 
erate double-well potential. As argued by Campbell et 
ai, breathers should form when the kinks have enough 
time to lose energy through their interaction, adjusting 
to their new trapped state. It is thus expected that these 
bound states form for small relative velocities for the 
kink-antikink pair, although the dependence on the ve- 
locity is far from trivial 1^. Once the breather forms, 
it will remain in its oscillatory state for a remarkably 
long time, with minimal emission of radiation. We are 
not aware of a detailed study of the lifetimes of these 
configurations for obvious reasons: it would take a huge 
amount of integration time in order to follow the evolu- 



tion of these bound states until their demise. The non- 
linear nature of the problem also has precluded (at least 
so far) analytical estimates for breather lifetimes. 

In the mid-1970's, 3-dimensional, breather-like config- 
urations were discovered in the context of 0^ models by 
Bogolubsky and Makhankov These authors found 
that, for certain initial conditions, spherically-symmetric 
bubbles settled into a long-lived configuration which they 
called pulsons. In 1994, one of us rediscovered these 
configurations while studying the collapse of subcritical 
bubbles in the context of degenerate and non-degenerate 
(p* models Since the typical behavior characteriz- 
ing these configurations is the high-frequency oscillation 
about the global minimum of the model, the name os- 
cillon was chosen instead. A further detailed analysis 
revealed more of the remarkable behavior of these con- 
figurations, such as the dependence of lifetime on ini- 
tial radius, the minimal radius for the bubbles to settle 
into oscillonic behavior, and the mechanism for their fi- 
nal demise . Recent work has pointed out the possible 
relevance of these structures for resonant hadronic states 

These configurations are far from being constrained 
to relativistic nonlinear field theories. Remarkably sim- 
ilar behavior has been found in experiments involving 
grains (or "sand") placed on a plate undergoing sinu- 
soidal vibrations 0. These PCFCs were independently 
named oscillons, and their discovery has triggered a host 
of theoretical work attempting to model the experimen- 
tal results. These include molecular dynamics simula- 
tions , semi-continuum theories 0, Ginzburg-Landau 
models [IC], coupled- map models [¥l|, order-parameter 
models [12|, and other continuum models |13|. The dif- 
ficulty here is in obtaining macroscopic laws describing 
the motion of granular materials, which can exhibit both 
solid and fluid properties. Oscillonic behavior has also 
been found in studies of acoustic instabilities in stars 
iQ. Long-lived, spatially-extended oscillatory behavior 
appears to be bringing together research in traditionally 
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very distant fields of physics. 

Our goal in this paper is to present a method which 
is extremely useful in the study of PCFCs arising in the 
context of nonlinear wave models. The numerical chal- 
lenge arises due to the peculiar nature of these configu- 
rations; although they are well-localized, they do radiate 
some of their energy, which, for small lattices, will get 
reflected back, compromising the numerical results. Fur- 
thermore, a large amount of energy may be shed initially 
as the field (or fields) settles (settle) into a PCFC. The 
simplest approach is to set up a very large spatial lattice, 
large enough that the outward radiation will never hit 
the boundary within the integration time. Clearly, this 
approach (used in Q) is extremely inefhcient for very 
long-lived PCFCs. An improvement is to use dynami- 
cally expanding lattices, that is, lattices that grow ahead 
of the radiation ||]. This saves some time, but not much 
as the lifetime becomes fairly large. Clearly, for more de- 
tailed studies of these objects, a more efficient numerical 
method is badly needed, one that allows for an accurate 
study of very long-lived PCFCs with relatively small lat- 
tices. Although we will introduce the method within the 
context of d-dimensional spherically-symmetric models, 
it can be easily implemented in more complex situations. 

The paper is organized as follows. In the next section 
we present in detail the numerical integration routine we 
used, the 4th-order operator splitting method developed 
in Ref . . This is followed by a discussion of the adi- 
abatic damping method we propose for studying long- 
lived configurations in small lattices. In section 3 we use 
3-dimensional oscillons to test the accuracy of the adi- 
abatic damping method by comparing its performance 
with results obtained from dynamically expanding lat- 
tices. In section 4 we present an analytical estimate for 
the minimum radius for the onset of 2-dimensional oscil- 
lonic behavior. This estimate is then tested numerically 
in section 5, where we also display the result for the life- 
time of a few oscillons, which live at least 3 orders of 
magnitude longer than their 3-dimensional counterparts. 
We conclude in section 6 with a summary of our results 
and directions for future work. 



II. NUMERICAL TECHNIQUES 
A. 4th order Symplectic Method 

The Hamiltonian field equations for our theory are 

^ = 9^0+ (l_l)a^0_a^y(0) (1) 

and 

^ = 7T , (2) 



where d is the number of spatial dimensions, and V{4>) 
is the potential, taken to be a function of the field 0. To 



integrate these equations, we use a higher-order operator 
splitting method, which is symplectic for Hamiltonian 
systems. In brief, symplectic methods use the fact that 
Hamiltonian systems of equations can be written as 



(3) 



where the vector z ~ (71^,0^) (in our case, z = 1, so 
we drop the i's), and where {a, 5} is a Poisson bracket. 
Now, we define the operator Dh by Dhz = {z, H}. The 
equations become 

z^Dhz , (4) 

which, integrating for time At, has the formal solution 

zit + At) ^ e-^*"" z{t) . (5) 

Note that, in general, systems of first-order equations can 
also be written in the form (||). However, the phase-space 
behavior of non-symplectic systems is usually singular. 
In general, Dh is a sum of terms Dh — Di+D2+---+Dn. 
To some order in At, we can approximate the exponential 
of a sum of operators as a product of exponentials, where 
each exponential has one of the operators as its argument. 
We use the notation 



to represent 



and 



to represent 



(At) 



(AtY 



/ AtDjv gAt_D2gAtDi^ 



(6) 



(7) 



(8) 



(9) 



So, for example, the 2nd order method 



,AtD, At_D2 „AtDi, 



)0 



,At£)A 



„AtD2„AtD 



is represented by 



(At){Atf 



e^*^i) (10) 



(11) 



In this paper, we use an explicit fourth-order method for 
'splitting' the operator exp(Z)//A<). The method is [Esl 



{Atf{At){Atf{-2At){Atf{Atf 
{Atf{Atf{At){Atf{At){At) 
{At){At){~2Atf{At){Atf{At) . 



(12) 



We split the Hamiltonian H = + (Vi/))^ + V{(t)) into 
two parts. Hi = (Vcj))^ + and H2 = tt^. The action 
of the operator exp{AtD2) on z is 
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0={<^,F2} = {0,ff} 
— IT 

7T = (13) 

which corresponds to integrating equation ( p^ 

'' + At7r" , (14) 



where the superscript n indicates the timestep, while 
leaving the value of tt the same. Note that this inte- 
gration is exact for these equations. The action of the 
operator exp(Ai£)i) on z is 



(/) = 



id -I) 



(15) 



which corresponds to integrating equation (|16| 



7r"+i = tt" + At{d 



(d-l) 



dp^ - d^V") (16) 



(we use fourth-order spatial differences for the derivative 
terms), while leaving the value of the same. Higher- 
order methods combine the integration of these equations 
in such a manner as to cancel higher-order errors in the 
commutators [Di,D2]- 

The standard leapfrog method is a second-order sym- 
plectic method. To see this from our notation, we write 



AtD2 ^^AtDi 



(17) 



where Diir = {tt, H} and D2(l) — {0, H}. Over the course 
of a simulation, we have 

e'«''e'^"''e'^"'') ^„„,„, . (18) 

Therefore, by first putting the momenta at the half- 
timestep, the leapfrog method is to alternately swap the 
integration of the 0's and tt's, then finally, to correct the 
momenta by a half-timestep. Higher-order methods such 
as the fourth-order method above require both swapping 
and integration with negative timesteps. 

The expected error in and tt from our method is 
RNAt^, where i? is a coefficient of order 0(1) given by 
the particular method and N is the number of timesteps 
in the simulation ||l5|. To test the above fourth-order 
method, we run a simulation with N — 650000 timesteps, 
where At — 0.01 on an expanding grid. Our expected er- 
ror at the end of the simulation is of order 6.5 x lO^^i?. 
To measure the actual accuracy of the method, we calcu- 
late the change in energy SE/E over the duration of the 
simulation and find a value 5E/E = 4 x lO^"' over 6500 
time units, or 650000 timesteps, which is compatible with 
the expected error. 



B. Adiabatic Damping Method: Small Lattices for 
Long Simulations 

In order to decrease the size of the grid, we want to 
absorb the non-linear radiation propagating away from 
the oscillon. For a problem involving massless radiation, 
we could use absorbing boundary conditions. However, 
for the problem at hand, the radiation is nonlinear, with 
a different dispersion relation for different wavemodes. 
Therefore, we introduce a damping term at a sufficient 
distance away from the oscillon such that the evolution of 
the oscillon is not disturbed, while the outgoing radiation 
is absorbed. This gives us a new equation of motion 



' + -f{p)4)-'sj'^(t> + V'{<t>)^Q . 



(19) 



The introduction of the non-zero decay term gives us a 
new equation for tt 



n+l 



1 



7(p)At/2 Vl-7(p)At/2 

{d-l 



-At(& 



2 J." 



P 



(20) 



Here, 7 is a function of p. Wc set 7 = for all grid- 
points less than pq. For gridpoints greater than pq, we 
set 7 = Tf'{p — Po)^, where 77 is a small constant. For 
the successful implementation of the method, it is cru- 
cial that po be chosen far enough away from the PCFC 
so as to not interfere with its dynamics. We find that, 
for a typical linear scale characterizing the PCFC of Rq, 
a safe choice is po ^ 20i?o, although different situations 
may call for different choices of po- In figure 1, we present 
a diagram of the relevant scales for the partition of the 
spatial lattice. 




S V 

FIG. 1. A schematic diagram for the relevant spatial scales 
in the implementation of the ADM. Rq is the linear scale of 
the PCFC, Rs is the location of the shell where the energy 
of the PCFC is measured, and po is the starting point of the 
damping term. 

The introduction of the damping term reduces our 
method from a symplectic method to an operator split- 
ting method, since the phase space of the system no 
longer obeys Liouville's theorem. 
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The smallness of r] ensures the very slow increase of 
the effective damping with the radial direction. (Hence 
the name adiabatic damping method.) This method can 
be easily generalized for higher dimensional lattices of 
different geometry. The p dependence of 7 was chosen 
such that the first derivative was zero at po. This choice 
gave us better accuracy than, for instance, a constant or 
an exponential dependence on radius. 

We should remark that the adiabatic damping equa- 
tions are mixed hyperbolic-parabolic. It has been claimed 
that since all higher-order operator splitting meth- 
ods have operators with backwards evolution in time, 
operator-splitting methods of order greater than 2 are 
unstable for parabolic problems. Clearly, this is incor- 
rect. We were able to integrate our mixed hyperbolic- 
parabolic equations for extremely long times (over 10^ 
timesteps) and found no instability. It is also clear why 
these methods are stable. The total backwards time inte- 
gration in our method is 4A<, but the method integrates 
forward in time for 16Ai. Therefore, any instability aris- 
ing from exponential increase in mode amplitudes from 
backward-time integration will be resuppressed by the 
forward-time integration. 



III. TESTING THE ADIABATIC DAMPING 
METHOD WITH 3D OSCILLONS 



The Lagrangian for our field theory is 
L = vr /" \2r]'-'^~^^dr - ^cj)'^ ~ T/(<^) 



(21) 



where d is the number of spatial dimensions (for ns d — 2 
or d — 3 only), and wc use the degenerate double- well 
potential. 



F(0) = 



A 



'^vacJ 



(22) 



As shown in , oscillons can easily be found by set- 
ting the initial field configuration with a Gaussian or a 
Tanh profile. Here, we will use the Gaussian ansatz, 

(j)= {(j)c- (pvac)e + (j)vac , (23) 

where Rq is the core radius and (l)vac is the asymptotic 
vacuum value of the field. (j)c sets the offset of the field 
from the vacuum at the core, the central displacement 
from equilibrium. 

The evolution of the configuration can be divided into 
three stages. Initially, the field sheds enough energy to 
settle (or not, if the initial parameters i?o and (j)c are 
outside the allowed range for oscillonic behavior) into 
the oscillon configuration. The lifetime of the oscillon 
stage is sensitive to the choices of Rq and (j)c, although 
the energy of the configuration is not. This is, in fact. 



what justifies identifying these PCFCs as a single config- 
uration. We conjecture that the lifetime of the oscillon 
configuration can be traced to the perturbations induced 
by the different choices of initial parameters, which will 
tend to increase the amount of radiation being emitted. 
However, we so far have not been able to prove this. The 
final stage is the oscillon's demise. As shown in Ref. [Q, 
due to the small but steady radiation from the oscillon, 
at some point the maximum amplitude allowed (that is, 
when (j)c = 0) falls approximately below the infiection 
point of the potential, the motion becomes linear, and 
(f'c — > (f'vac exponentially fast. 

As one would expect, one must understand something 
about the evolution of the oscillon to correctly set po, 
since damping must not interfere with the evolution of 
the oscillon. As a rule of thumb, the better localized 
the configuration, the closer po can be to Rq. The fact 
that oscillons behave asymptotically as exp(— p) certainly 
helps. As remarked above, we found that po = 20i?o 
worked well. 

To make sure that the damping region we have intro- 
duced does not interfere with the evolution of the oscillon, 
we checked that the oscillon lifetime is not significantly 
different on grids with damping. We simulate oscillons 
on two grids. One simulation is performed on an expand- 
ing grid (Aa; = 0.2), the other on a 1024 point grid with 
physical size 204.8 (Aa; — 0.2), where the damping with 
constant t] = 0.005 begins a.t p — 20i?o- In figure 2, we 
have superimposed the graph of an expanding grid sim- 
ulation, which we take to be the exact result, and the 
approximate result, which comes from a simulation with 
a damping region for several choices of initial radii, Rq. 
(All physical quantities are quoted in dimensionless units. 
For 3-dimensional relativistic field theories, distance and 
time are given in units of A"^/^^"^,^^,.) The plots denote 
the time evolution of the total energy in a shell about the 
core (we have set our shell at Rg = 5i?o)- 

As is apparent from figure 2, oscillon lifetimes are es- 
sentially indistinguishable for the time spacings of the 
data output in these graphs. Evolution of oscillons gen- 
erated from initial configurations with larger radii agrees 
less well with the exact oscillons due to a greater inter- 
action between long-wavelength modes, generated as the 
oscillon decays at the end of its lifetime, and the reflecting 
boundary at the far end of the damped part of the grid. 
This can be cured by using a larger grid, and increasing 
Pq. Notwithstanding these considerations, the agreement 
in lifetimes is still extremely good for 1024 gridpoints, as 
is shown next. The reason that long- wavelength modes 
introduce difficulties with this method is that they are 
damped less efficiently than short-wavelength modes (the 
damping is roughly proportional to exp{—k'^t)) . 
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IV. ANALYTICAL ESTIMATE FOR ONSET OF 
2-DIMENSIONAL OSCILLON 
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FIG. 2. Shell energy vs. Time: The evolution of the energy 
in a shell of radius Rs = 5Ro about the oscillon core as a func- 
tion of time. Two graphs are superimposed, the exact solution 
and the solution using a small 1024 point grid with damping. 
From left to right, we see the curves for Ro — 2.5, 2.6, and 2.7. 



In figure 3, we plot the logarithm (base ten) of the 
absolute value of AEg/Eg as a function of time, where 
Eg is the energy in a shell of radius 5i?o about the origin. 
We chose Rq — 2.7 in this example. AEg is the difference 
between the energies of the expanding grid simulation 
and the damped simulation. The shell energies agree 
with each other to better than one hundredth of a percent 
until the oscillon decays. 
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FIG. 3. AEs/Es vs. Time: The relative difference between 
the expanding grid and damped grid solutions as a function 
of time for _Ro = 2.7. Note that through all of the lifetime of 
the oscillon the solutions agree to better than one hundredth 
of a percent. 



In order to estimate the approximate minimum size 
at which oscillon-like PCFCs appear in a 2-dimensional 
nonlinear field theory, we begin with the ansatz 



{r,t) = A{t)e ^ + . 



(24) 



where A{t) = 4>c{t) — 4>vac- The Lagrangian for our field 
theory is given in ( pT| ) with d — 2, and the potential is 
given in ( p^ as in the case for 3d oscillons. This ansatz 
restricts the time-dependence to the amplitude at the 
core of the configuration, taking the radius to be a con- 
stant parameter. Despite its simplicity, reducing a field 
theory to a model with one degree of freedom, it was 
used successfully for determining the minimum radius of 
3d oscillons in Ref. (5). 

Substituting the ansatz (|2^ ) into the Lagrangian and 
integrating over the radial dimension gives an effective 
Lagrangian for the single degree of freedom A{t), 



L = 2tt 



Rn 



Rn 



Rn 



8 4 32 6 4 



The equation of motion is 



Rl 



A + 2A- 



1 



-A'^ -2A^ = 



(25) 



(26) 



We then assume that we can separate A{t) = Ao{t) + 
6A{t) and investigate harmonic perturbations about the 
background (assumed stable) solution Ao{t). This leads 
to a frequency response of 



u^Ao, R„) = 2+-^ + -Al- 4Ao 
Rq ^ 



(27) 



If we now minimize the frequency as a function of Aq, we 
find 



~3 



(28) 



We expect an instability to develop for oj'^ < 0, giving a 
bifurcation at Rq = \/3: oscillons can only exist in the 
presence of these instabilities, which ensure the (tempo- 
rary) survival of the nonlinear regime. This is reminis- 
cent of the well-known spinodal instability in Ginzburg- 
Landau systems, where the growth of instabilities occurs 
as the system probes the concave part of the potential 
(or free-energy density) ||l^. In short, although the per- 
turbed solution is exponentially increasing for this linear 
analysis, we expect nonlinear terms to stabilize it, at least 
temporarily, leading to oscillon solutions. As we will see 
below, this prediction is quite accurate. 
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V. RESULTS FOR TWO-DIMENSIONAL 
OSCILLONS 

We can now use the adiabatic damping method to 
study 2d oscillons, which preliminary studies (with large 
lattices) have shown to be remarkably long-lived. (Pre- 
vious simulations were stopped as lifetimes went over 
T 10'* time units.) In figure ^, we plot the energy 
in a shell of radius Rs — 5i?o about the origin as a func- 
tion of time. Starting with Gaussian initial configura- 
tions, we searched for oscillons around Rq — -s/S, the 
expected onset value of the oscillon solution obtained 
above, and found the bifurcation at the slightly lower 
value of Rq ~ 1.71, as indicated in the figure. 
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FIG. 4. Shell energy vs. time: The evolution of the energy 
in a shell of radius Rs = 5Ro about the oscillon core as a 
function of time. Initial conditions were, from bottom to top, 
Ro = 1.69, 1.71, 1.72, 1.73, and 1.75. 

Above this value for the initial radius, our results ex- 
hibit the same behavior as for the 3d oscillons, with a few 
new and surprising properties, which are easily identified 
in figure First, we could not find a maximum radius 
above which the initial configurations do not become os- 
cillons. (We searched all the way to Rq = 30.) We recall 
that for 3d oscillons obtained from Gaussian initial con- 
figurations, the maximum radius was Rq ~ 4.2 [||. We 
suspect this peculiar behavior to be a consequence of the 
fact that, in 2d, Gaussian configurations have constant 
surface energy densities, that is, no radial dependence. 
[For 2d Gaussians, the energy goes a.s E = A + BR^, 
where A and B depend on the potential.] For initial con- 
figurations with Tanh profiles, we found oscillons within 
the interval 1.5 ^ Rq ^ 4.5, as indicated in figure ^. We 
note that Tanh configurations with small radii are fairly 
well approximated by Gaussian configurations, that is, 
have small surface energy terms. 
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FIG. 5. Shell energy vs. time: Examples of several 2d os- 
cillons. The radii for the initial Gaussian profiles are denoted 
by each respective curve. 

Second, the energy plateau is not as clearly defined 
as in the 3d case. For example, a configuration with 
initial radius Rq = 10 evolves into a lower energy oscil- 
lon than the other cases in the figure. Also, the oscillon 
energy for the marginal case Rq = 1.73 is lower than 
the rest. This suggests that there may exist several "ex- 
cited", or metastable, oscillonic states in 2-dimensional 
models, something worth pursuing. 




1.5 

- 5. 




t I I I I I I I ^ ^ I I ^ I I I; 

1000 2000 3000 

Time 

FIG. 6. Shell energy vs. time: Examples of several 2d 
oscillons from a Tanh initial configuration. The initial radii 
are denoted by each respective curve, with Ro = 2, and 3 
omitted for clarity. 

Third, we also could not find a finite lifetime for the 
oscillons. Using the adiabatic damping method, we fol- 
lowed the evolution of the oscillons for over 10^ time units 
without observing their demise. [With our time steps this 
implies over 2 x 10^ time-iterations.] This strongly sug- 
gests that these may be time-dependent stable solutions. 
We are currently pursuing this issue in more detail, us- 
ing both analytical and numerical methods. We note 
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that in theories with conserved particle number, simple 
time-dependent solutions known as non-topological soli- 
tons have been found However, in the case at hand, 
the conserved quantity giving rise to the (possible) sta- 
bility of the configuration is not immediately obvious. 

VI. SUMMARY AND OUTLOOK 
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We have presented a method designed for the numeri- 
cal investigation of long-lived field configurations such as 
breathers, oscillons and other spatially-extended persis- 
tent coherent field configurations (PCFCs). The method 
uses an adiabatically increasing damping term in the 
equation of motion, placed safely away from the PCFC 
so as not to interfere with its dynamics. We argued that, 
with this method, it is possible to follow the dynamical 
evolution of these objects for extremely long times, allow- 
ing us to obtain accurate results with very small lattices. 
Using 3-dimensional oscillons as an example, we showed 
that the method allows for accuracies of a fraction of a 
percent in the measurement of physical quantities, such 
as the energy and lifetime, typically of the order of a 
part in 10^. We then applied this approach to inves- 
tigate 2-dimensional oscillons. After obtaining the mini- 
mum radius that allows for their existence, we discovered 
that lifetimes can exceed 10^ time units. Although not 
yet proven, it is possible that these time-dependent field 
configurations are absolutely stable. 

There are quite a few obvious avenues for future work, 
in addition to the ones mentioned above. It would be 
straightforward to apply the method to more complex 
geometries, for example searching for tube-like oscillons, 
or the study of interactions between oscillons, similar to 
the I-dimensional kink-antikink scattering studies. As 
mentioned elsewhere, these configurations may be ther- 
mally nucleated during phase transitions, leading to im- 
portant corrections to decay rates and completion times 
Given the longevity of 2d oscillons, this may be par- 
ticularly relevant to systems in the 2-dimensional Ising 
universality class. It would also be interesting to study 
"excited oscillons" , that is, non spherically-symmetric 
configurations which can be expanded into a series of har- 
monics. Do these configurations decay into their ground 
state (a normal, £ — oscillon) or are they completely 
unstable? Hopefully, future research will establish the 
connection between the various oscillons described in the 
literature as obeying some simple general principles. 
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